
%% initialize
clear all;
close all;

rng(328120)

Nsim = 500;

loadstarting = 0;
loadestimates = 0;

%% read in data

% choices

rawdata = csvread('../../01_data/02_basedata/estimationdata/All_T_J/All_T_J_choices_buff_0_0_lf_principalExtra_C.csv',1,0);

IJjtidx         = rawdata(:,1);
IJsidx          = rawdata(:,2);
IJpositive      = rawdata(:,3);
IJVAmadev       = rawdata(:,4);
IJVAmean        = rawdata(:,6);

IJX             = [IJVAmadev+IJVAmean,rawdata(:,5)+rawdata(:,7),rawdata(:,8:end)];

clearvars rawdata;

rawdata = csvread('../../01_data/02_basedata/estimationdata/All_T_J/All_T_J_choices_buff_0_0_lf_principalExtra_cf_C.csv',1,0);

IJitidx         = rawdata(:,1);
IJjtidx_teach   = rawdata(:,2);
IJjtidx_princ   = rawdata(:,3);
IJapp_year_cf   = rawdata(:,4);
IJVAmadev_cf    = rawdata(:,5);
IJVAmean_cf     = rawdata(:,7);

IJX_princ       = [IJVAmadev_cf+IJVAmean_cf,rawdata(:,6)+rawdata(:,8),rawdata(:,9:end)];


clearvars rawdata;




%% prepare data for estimation

ds = max(IJsidx);
d2 = max(IJjtidx);



%% draw random effects


Xdraws = zeros(size(IJjtidx,1),1,Nsim);

drawss = normrnd(0,1,[d2,Nsim]);

Xdraws = zeros(size(IJjtidx,1),1,Nsim);
Xdraws(:,1,:) = drawss(IJjtidx,:);

%% estimate mixed logit


y_logistic = IJpositive;


Xmat = [ones(length(IJpositive),1),IJX];

f =  @(x)func_mixed_logit(x,y_logistic,Xmat,Xdraws,IJjtidx,0);


    betastart = zeros(size(Xmat,2),1);
    sigmastart  = ones(size(Xdraws,2),1);
    
    betastart(1) = -1;
 
    thetastart = [betastart;sigmastart];




optionsgrad = optimset('MaxFunEvals',10000000,'MaxIter',10000000,'TolX',1E-5,...
    'Display','Iter','TolFun',1E-7,'Algorithm','trust-region','GradObj','on');


    tic
thetahat = fminunc(f,thetastart,optionsgrad);
toc
f =  @(x)func_mixed_logit(x,y_logistic,Xmat,Xdraws,IJjtidx,1);

[obj,grad,se] =  f(thetahat);

 save '../../01_data/05_temp/mixed_logit_principal_totVA_C'
 

